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We analyze a specific class of random systems that are driven by a symmetric Levy stable noise, 
where Langevin representation is absent. In view of the Levy noise sensitivity to environmental 
inhomogeneities, the pertinent random motion asymptotically sets down at the Boltzmann-type 
equilibrium, represented by a probability density function (pdf) p*{x) ~ exp[— $(x)]. Here, we 
infer pdf p(x,t) based on numerical path-wise simulation of the underlying jump-type process. A 
priori given data are jump transition rates entering the master equation for p(x,t) and its target 
pdf p*(x). To simulate the above processes, we construct a suitable modification of the Gillespie 
algorithm, originally invented in the chemical kinetics context. We exemplified our algorithm simu- 
lating different jump-type processes and discuss the dynamics of real physical systems where it can 
be useful. 

I. INTRODUCTION 

Despite many attempts to pin down the essential features of dynamics and relaxation in random systems, the 
problem is still far from its complete solution. It turns out that Levy flight models are adequate for the description 
of different random systems ranging from the motion of defects in disordered solids to the dynamics of assets in 
stock markets, see, e.g. They are especially useful to model the random systems on the semi-phenomenological, 
mesoscopic level, when the (often unknown) details of their microscopic random behavior are substituted by a properly 
tailored (e.g. based on experimental data) (Gaussian or Levy) noise. Paradoxically, in disordered solids, the noise can 
promote order and organization, switching them between different equilibrium states. The latter situation emerges 
in disordered ferroelectrics, where the fluctuations of order parameter (spontaneous polarization) give rise to self- 
localization of charge carrier, generating a nuctuon, an analog of well-known polaron in disordered substance 0]. 
These nuctuons make a substantial contribution into conductivity and optical properties of disordered dielectrics. 

Many random processes admit a description based on stochastic differential equations. In such case there is 
a routine passage procedure from microscopic random variables to macroscopic (statistical ensemble) data. The 
latter are encoded in the time evolution of an associated probability density function (pdf) which is a solution of a 
deterministic transport equation. A paradigmatic example is the so-called Langevin modeling of diffusion-type and 
jump-type processes. The presumed microscopic model of the dynamics in external force fields is provided by the 
Langevin (stochastic) equation whose direct consequence is the Fokker-Planck equation, Q and 0]. We note that 
in case of jump-type processes the familiar Laplacian (Wiener noise generator) needs to be replaced by a suitable 
pseudo-differential operator (fractional Laplacian, in case of a symmetric Levy-stable noise). 

We pay a particular attention to jump-type processes which are omnipresent in Nature (see @ and references 
therein). Their characterization is primarily provided by jump transition rates between different states of the system 
under consideration. In the present paper, our major focus is on a specific class of random systems which are plainly 
incompatible with a straightforward Langevin modeling of jump-type processes and, as such, are seldom addressed in 
the literature. 

To this end we depart from the concept, coined in an isolated publication Q, of Levy flights-driven models of disorder 
that, while at equilibrium, do obey detailed balance. The corresponding research line has been effectively initiated in 
Refs. It has next been expanded in various directions, with a special emphasis on so-called Levy-Schrodinger 

semigroup reformulation of the original probability density function (pdf) dynamics, [1, [TH, [l9[ and [l(|-[lH, c.f. also 
jl7Hl9| . We note in passing that the familiar Fokker-Planck equation can be equally well reformulated in terms of the 
Schrodinger semigroup and this property is universally valid in the standard theory of Brownian motion, 0, [TtJ • Its 
generalization to Levy flights is neither immediate nor obvious. It is often considered in the prohibitive vein following 

HUSH- 

In fact, in relation to Levy flights, a novel fractional generalization of the Fokker-Planck equation has been intro- 
duced in Refs. Q- @ to handle systems that are randomized by symmetric Levy-stable drivers. In this case, contrary 
to the popular lore about properties of (Langevin- based) Levy processes (c.f. Refs. [2(|- HH and HH), the pertinent 
random systems are allowed to relax to (thermal) equilibrium states of a standard Boltzmann-Gibbs form. 

The underlying jump- type processes, in the stationary (equilibrium) regime, respect the principle of detailed balance 
by construction [la, Hal- Their distinctive feature, if compared with the standard Langevin modeling of Levy flights, 
is that they have a built-in response not to external forces but rather to external force potentials. These potentials are 
interpreted to form confining "potential landscapes" that are specific to the environment. Levy jump-type processes 
appear to be particularly sensitive to environmental inhomogeneities, fA [Til]. 
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We note in passing that confining Levy nights resolve themselves to the problem of truncated distributions, which 
can be addressed in many ways, beside a simple cut-off, by considering fast falling tails. The random walk theory 
predominantly employs to the master equation as the major tool to quantify random motion. Inhomogcncous tran- 
sition rates have been used in the literature before and the emergence of non-trivial target distributions has been 
demonstrated, c.f. [24j . 

Levy flights are pure jump (jump-type) processes. Therefore, it seems useful to indicate that various model real- 
izations of standard jump processes (jump size is bounded from below and above) can be thermalizcd by means of a 
specific scenario of an energy exchange with the thermostat. It is based on the principle of detailed balance. We have 
discussed this issue in some detail before along with an extension of this conceptual framework to Levy-stable 
processes. Not to reproduce easily available arguments of past publications, we shall be very rudimentary in our 
motivations, see however [jq . 

We quantify a probability density evolution, compatible with a jump-type process on R (this limitation may in 
principle be lifted in favor of R n ), in terms of the master equation: 



d t p(x,t) = j [w <l> {x\y)p{y,t)-w lt> {y\x)p{x,t)]dy, (1) 

El<|x— J/|<£2 

where £\ and £2 are, respectively, the lower and upper bounds of jump size and 

exp[($(y) - *(a0)/2] 



w<t>( x \y) = c n- 



- ?y| 1 +/ i 



\x - y\ 



c = L(l + p) sin(7r/i/2) 

^ TT 

is the jump transition rate from y to x. We stress that w^,(x\y) is a non-symmetric function of x and y. 

An implicit Boltzmann-type weighting involves a square root of a target pdf p*(x) ~ exp[— $(2;)] and accounts 
for the a priori prescribed "potential landscape" $(x) whose confining features affect the jump-type process. What 
matters is a relative impact of a confinement strength of $(x) (level of attraction, see Ref. |9|) upon jumps of the 
size \x — y\, both at the point of origin y and that of destination x. In principle, $>(x) may be an arbitrary function 
that secures a i 1 (fl) normalization of exp(— $(ar)). In this case, the resultant pdf p* is a stationary solution of the 
transport equation (p} with unbounded jump length, e.g. E\ — > and £2 — > 00. 

We note that the presence of lower and upper bounds of the jump size £1,2, that are necessary for an implementation 
of numerical algorithms, enforce a truncation of the jump-type process (without any cutoffs) to a standard jump 
process. The transition rates of the latter, however, are ruled by Levy measures of symmetric Levy stable noises 
with p S (0,2). A lower bound for the jump size is usually removed while evaluating the corresponding integrals in 
the sense of their Cauchy principal values. An upper bound is less innocent and its effects need to be controlled by 
long tailed pdfs which stands for a distinctive feature of Levy flights, see a discussion of Levy stable limits of step 
processes in Ref. flCfl. There is also pertinent discussion of a long time behavior of (unconfined, e.g. free) truncated 
Levy flights in RefHS- 

In contrast to procedures based on the Langevin modeling of Levy flights in external force fields, [1, HH, H^, there is 
no known path- wise approach underlying the transport equation ([1]). With no direct access to sample trajectories of 
the stochastic process in question, a method must be devised to generate random paths directly from jump transition 
rates ([5]). The additional requirement here is that we set a priori a "potential landscape" $(x) for a chosen jump-type 
(symmetric Levy stable) noise driver. 

The outline of the paper is as follows. First we describe our modification of the Gillespie algorithm which entails a 
numerical generation of random paths for the dynamics determined by Eqs (flj and @ . Next the statistics of random 
paths is addressed and various accumulated data are analyzed with a focus on inherent compatibility issues. 

We analyze generic (Cauchy, quadratic Cauchy) and non-generic (Gaussian and locally periodic) examples of target 
pdfs for the jumping dynamics. Random paths are generated in conjunction with representative Levy stable drivers, 
like e.g. those indexed by p — 1/2, 1, 3/2. Their qualitative typicality is emphasized. 

Statistical data, acquired from our modification of Gillespie algorithm, have been employed to generate the dynam- 
ical patterns of behavior p(x,t) — > p*(x). Effectively, that entails a path-wise reconstruction of the solution p(x,t) 
of the master equation |T]), ([2"]l. or - in case this solution it is available prior to trajectory simulations - to verify a 
compatibility of the transport (master) equation (flj , © and its underlying path- wise representation. All that comes 
from the predefined knowledge of the target pdf and non-symmetric (biased) jump transition rates. 
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II. MODIFIED GILLESPIE ALGORITHM. 

Here we adopt [29| (and properly adjust to handle Levy flights) basic tenets of so-called Gillespie's algorithm 
[26J. Originally, this algorithm had been devised to simulate random properties of coupled chemical reactions. 
The advantage of the algorithm is that it permits to generate random trajectories of the corresponding stochastic 
process directly from its (jump) transition rates, with no need for any stochastic differential equation and/or its explicit 
solution. We emphasize that this feature of Gillespie's algorithm is vitally important, since Langevin modeling is not 
operational in our framework. 

We rewrite Eq. ([T]) in the form [x — y = z) 



d t p(x, t) 



w^{x\z + x)p(z + x, t) w c j 1 (z + x\x)p(x, t) 



dz. 



(3) 



e 1 < I z I < e 2 



To construct a reliable path generating algorithm consistent with Eq. ([3]) we first note that chemical reaction channels 
in the original Gillespie's algorithm may be re-interpreted as jumps from one spatial point to another, like transition 
channels in the spatial jump process. An obvious provision is that the set of possible chemical reaction channels is 
discrete, while we are interested in continuous set of all admissible jumps from a chosen point of origin xq. With a 
genuine computer simulation in mind, we must respect standard numerical assistance limitations. Surely we cannot 
admit all conceivable jump sizes. As well, the number of destination points, even if potentially enormous, must remain 
finite for any fixed point of origin. 

Our modified version of the Gillespie's algorithm, appropriate for handling of spatial jumps is as follows [30j |: 

1. Set time t = and the point of origin x = xq. 

2. Create the set of all admissible jumps from xq to xq + z that is compatible with the transition rate w^(z + xq\xo)- 



3. Evaluate 



Wx{x ) 
W 2 (x ) 



w<j,(z + x \x )dz, 

2 

w<j,(z + x \x )dz 



(4) 



and W(xo) = W 1 (x ) + W 2 (x ). 

4. Using a random number generator draw p S [0, 1] from a uniform distribution. 

5. Using above p and identities 



/ Wj,(z + x \x )dz = pW(x ), 



p< WM/Wixo); 



Wi(x ) + J w^z + x \x )dz = pW(x ), p > W 1 (x a )/W{x ), 



(5) 



find b corresponding to the "transition channel" xq — > b. 

6. Draw a new number q G (0, 1) from a uniform distribution. 

7. Reset time label t = t + At where At = — \nq/W(xo)- 

8. Reset xq to a new value xo + b. 

9. Return to the step 2 and repeat the procedure anew. 



III. STATISTICS OF RANDOM PATHS FOR DIFFERENT JUMP- TYPE PROCESSES 



In this section, we select the jump- type processes, which best demonstrate the capabilities of our algorithm of random 
paths generation. Once suitable path ensemble data are collected, we ca n reconstruct the pdf p{x, t) dynamics and next 
verify whether statistical (ensemble) features of generated random trajectories are compatible with the asymptotic 
a priori imposed upon solutions of the master equation (|T|). That refers to a control of an asymptotic behavior 
p(x, t) — > (x) when t — > 00 . 



4 




FIG. 1: Gaussian target: Time evolution of the pdf p(x,t) second moment for 75 000 trajectories. Inset visualize the oscillations 
smoothing in the asymptotic regime for 10 < t < 15; figures near curves correspond to p values. 




FIG. 2: Gaussian target: Time evolution of p(x,t) inferred from 75 000 trajectories: p = 0.5 (left panel), p — 1 (middle panel) 
and p = 1.5 (right panel). All trajectories originate from x — 0, i.e. refer to the 8(x)-type initial probability distribution. 



A. Harmonic confinement 



Let us first consider an asymptotic invariant (target) pdf in the Gaussian form: 



p*(x) = ^=e-*\ (6) 

\/7T 



The corresponding /^-family of transition rates reads 

£ -z 2 /2+xz 

W4> (z + x\x) = C M ^ . (7) 

The C-codes for trajectory generating algorithm, |30| . were employed to get the trajectory statistics data for Levy 
drivers with p, = 0.5, 1, 1.5. 
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The results of our numerical simulations are reported in Figs. 1 and 2. We note, that in Fig. 1 the second 
moment oscillates near its equilibrium value 1/2. A numerical convergence to < X 2 >=l/2 is consistent with an 
analytic equilibrium value of the second moment of the chosen p*(x) El The rate of this convergence is higher for 
larger p € (0, 2). Clearly, for small p the big jumps are frequent which enlarges the inferred time intervals At in the 
Gillespie's algorithm. Thus, the relaxation to equilibrium is slow. It gets faster for larger p, when big jumps are rare 
and time intervals At are generically very small. 

Fig. 2 displays a probability density evolution, inferred from the ensemble statistics of 75000 trajectories. All of 
them have started form the same point x = 0. Although the data fidelity grows with the number of contributing 
paths, we have not found significant qualitative differences to justify a presentation of data for 100 000, 200 000, 250 
000 and more trajectories. The relaxation time rate dependence on p is clearly visible as well. It suffices to analyze 
differences between three curves for t = 0.2 and/or t = 1. We observe a conspicuous lowering of their maxima with 
the growth of p (take care of different scales on the vertical axes on Fig. 2 panels). The simulated pdfs at t = 10 
are practically indistinguishable from an exact analytical asymptotic pdf El The convergence of p{x,t) towards p*(x) 
appears to be relatively fast irrespective of the chosen //-driver. 

Although our reasoning is definitely path- wise and all data have been extracted from simulated trajectory ensembles, 
our focus was on the inferred dynamics of p(x,t). We do not reproduce figures visualizing generic sample paths. 
However, we indicate their (paths) distinctive qualitative typicalities, if one would compare e.g. motion scenarios 
corresponding to stability indices p = 1/2, p = 1 and p = 3/2, The structural impact (probability of occurrence) 
of larger against smaller jumps, even on the visual level, conforms with standard simulations of Levy stable sample 
paths (with no forces or potentials involved), c.f. [3lj . 



B. Logarithmic confinement 



1. Quadratic Cauchy target 



Let us consider a long-tailed asymptotic pdf which is a special a — 2 case of the one-parameter a-family of 
equilibrium (Boltzmann-type) states, associated with a logarithmic potential = aln(l + a; 2 ), a > 1/2, see 

\MM ■■ 

2 1 

7T (1 + X i Y 

The transition rate[2] w$(z + x\x) for any p € (0, 2) takes the form 

W ^ Z + xlx)= \z C \^l + ('z X +xy (9) 

Simulation results for this case are displayed in Figs. 3 and 4. If we compare Fig. 3 with Fig. 1 we see the 
existence of small oscillations in the asymptotic regime about the value 1/2. Those from Fig.l were relatively small, 
while those on Fig. 3 are more noticeable. This is related to much slower decay of transition rates ([9]) (determined 
by slow-decaying asymptotic pdf flSJ)), as compared to those for Gaussian case J?]). 

The second moment of the present p*(x), ([5)1, equals 1 and the convergence towards this value is clearly seen in 
Fig. 3. This convergence is much slower than in the Gaussian (harmonic confinement) case which is not a surprise: 
dZJ) and (|9]) indicate that the present rate of convergence should be logarithmically slower. Fig. 5, quite alike Fig. 2, 
convincingly demonstrates a convergence of p{x, t) to the asymptotic p*(x). For definitely large times around t = 400, 
p(x, t) and p*{x) become practically indistinguishable. Similarly to the Gaussian case, the rate of convergence becomes 
larger with the growth of p £ (0, 2). 

2. Cauchy target 
Now we consider an asymptotic pdf of the form : 



P*( x ) = -TT-j- ( 10 ) 

TT 1 + X Z 



In this case, the transition rate from x to x + z reads 



/ 1 + x 2 



w ,(z + x\x) = w ^J TTI —y . (11) 
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FIG. 3: Quadratic Cauchy target: Time evolution of the pdf p(x,t) second moment for 200 000 trajectories. 
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FIG. 4: Quadratic Cauchy target: Time evolution of p(x,t) inferred from 200000 trajectories for n = 0.5 (left panel), fi — 1 
(middle panel) and p: — 1.5 (right panel). All trajectories are started from x — 0. Note scale differences on vertical axes. 



We consider Cauchy driver corresponding to /i = 1. In Fig. 5, we report the time evolution of the statistically inferred 
p(x, t) for different time instants. An approach to the asymptotic pdf ((TO)) is clearly seen, together with a convergence 
of a half- width to its asymptotic value 1. It can be shown that the same convergence pattern is valid for cumulative 
probability distribution which approaches the asymptotic function F(x) = h + arctan x , 



C. Locally periodic confinement 



To set firm grounds for future research it is instructive to study our model for more complicated forms of confining 
potentials. In view of their physical relevance, it is appealing to address an issue of confining (trapping) environments 
with a periodic spatial structure. This problem may be relevant to the random motion of defects in doped semicon- 
ductors and to phenomena like hopping conductivity. Here, we encounter a major difficulty with intcgrability of the 
Boltzmann-type weighting function exp(— $) as for periodic $ the corresponding integral is divergent. Periodicity 
and integrability can here be reconciled by means of locally periodic potentials that take a definite confining form 



FIG. 5: Cauchy target: Time evolution of pdf p(x,t) for j-i — 1 and 200000 trajectories, starting from x — 0. 



(harmonic or polynomial) for larger values of x G R. Let us consider the following asymptotic pdf 



1 e -sin 2 (2^) I I < 2 . 



C ' I I — ' 



P*{x) = 



ie-(- 2 - 4 ), |z| > 2, 



(12) 



where C = 3.032818 is a normalization constant 



The transition rate from x to x + z reads 



w ( /,(z + x\x) 




exp [{4>{x) - <p{x + z))/2 



(13) 



where the potential <f> has the following locally confining form 



4>{x) = 




(14) 



Time evolution of the inferred pdf p(x,t) is reported in Fig. 6 for \x = 1. All sample trajectories were started 
form x = which corresponds to the <5(cc)-type initial distribution. The probability density spreads out with time 
in conformity with the trapping (confining) properties of the locally periodic enclosure (environment or "potential 
landscape"). For large running times t=400 the trajectory statistics produces data that are indistinguishable from 
those for the asymptotic pdf. We have checked that beginning from about 100 000 trajectories, further accumulation 
of the trajectories number like e.g. 200 000 (displayed) and 300 000 (not displayed) for the data statistics is inessential. 
In such cases the curves are almost the same, we merely improve a fidelity of the statistics. 



If a random process does not admit the description in terms of a stochastic differential equation (e.g. Langevin 
modeling), its direct numerical simulation becomes impossible by means of existing popular algorithms. In the 
present paper, for the first time in the literature, we propose a working method to generate stochastic trajectories 
(sample paths) of a random jump-type process without resorting to any explicit (or numerical) solution of a stochastic 
differential equation. 

To this end we have modified the Gillespie algorithm (25l . [2r| , normally devised for sample paths generation, if 
transition rates occur between a finite number of states of a system. The essence of our modification is that we take into 



IV. CONCLUSIONS. 
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FIG. 6: Time evolution of p(x,t) inferred from 200000 trajectories at p = 1. The data for 100000 and 300000 trajectories (not 
displayed) do not show qualitative differences. 



account the continuum of possible transition rates, thereby changing the finite sums in the original Gillespie algorithm 
into integrals. The corresponding procedure for stochastic trajectories generation has been changed accordingly. 

In other words, here we "extract" the background sample paths of a jump process, whose pdf p(x, t) is to obey the 
generalized Fokkcr-Planck dynamicsHJEJ • We emphasize once more here, that we have focused our attention on those 
background jump-type processes that cannot be modeled by any stochastic differential equation of the Langevin type. 
However, our ultimate result was a reconstruction of the pdf dynamics from the path-wise data. Thus, our modfication 
of the Gillespie's algorithm serves two purposes: (i) getting access to path-wise data and (ii) reconstructing the pdf 
p{x, t) dynamics, compatible with the master equation in question, without resorting to other methods of its solution 

Although heavy-tailed Levy stable drivers were involved in the present considerations, we have clearly confirmed 
that an enormous variety of stationary target distributions is dynamically accessible in each particular p € (0, 2) case. 
That comprises not only a standard Gaussian pdf, discussed usually in relation to the Brownian motion like Wiener 
process. We can reach Cauchy pdf in terms of any p ^ 1 driver as well, provided a steering environment is properly 
devised. In turn, the Cauchy driver in a proper environment may lead to an asymptotic pdf with a finite number of 
moments, the Gaussian case included (jlal)- 

An example of the locally periodic environment has been considered as a toy model for more realistic physical 
systems. We mention possible generalizations of our method to the Brownian motor concept (see, e.g., Ref. [32| for 
recent review) to include a non-Gaussian jumping component. In those systems it is the properly tailored periodic 
" potential landscape" which enforces a conversion of a homogenous stochastic process (Brownian motion for reference) 
into the directed motion of particles at nanometer scales. It is closely related to the problem of so-called sorting in 
periodic potentials [33| . Other problems to be addresses concern ultracold atoms in optical lattices subject to random 
potentials |34| . which are promising not only from purely scientific point of view, but also for many technological 
applications. 

We note that the theoretical description of the above mentioned topics relies essentially on the Langevin-like equation 
input. Our approach offers prospects for generalizations, where systems with non-Langevin response to external 
potentials may come into consideration, along with more traditional ones. What we actually need to implement our 
version of Gillespie's algorithm is the knowledge of jump transition rates of those random systems only. 

Our preliminary work (in progress) shows that an extension of our algorithm to higher dimensions is operational. 
In particular, the planar case is worth further analysis, possibly with more realistic "potential landscapes". The direct 
numerical modeling in three-dimensional case can also be important, for instance, in the investigations of the dynamics 
of amorphous materials. Theoretical studies of their physical properties have a longer history, starting in the eighties 
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and even earlier [35|,[36j]. These materials may comprise conventional glasses (see, e.g., [13), so-called spin [HI,!!?} and 
orientational (for example, dipole) glasses [36[ as well as virtually any disordered solid like doped semiconductors. In 
these substances, the random motion of defects can be regarded as real jump-type process to be modeled numerically. 
Moreover, if defects have only rotational degrees of freedom (like spins or electric/elastic dipoles), their dynamics 
can be well represented by some effective jump-type process. Such representation can be complementary to existing 
Monte-Carlo modeling techniques which for disordered systems are very computer intensive. For instance, we hope to 
tackle the issue of interplay between exponential and long-time (logarithmic or stretchcd-cxponcntial) relaxation in 
these substances. Also, the 3D modification of our algorithm can be used in the simulation of switching processes in 
disordered ferroelectrics [38| . Similar fictitious jump- type process can be assigned to the dynamics of inhomogeneous 
broadening of resonant lines [39l ]. This broadening occurs in condensed matter and/or biological species, in a number 
of spectroscopic manifestations like the electron paramagnetic resonance, nuclear magnetic resonance, optical and 
neutron scattering methods. It arises due to random electric and magnetic fields, strains and other perturbations 
from defects in a substance containing the centers whose resonant transitions between energy levels are studied. Here 
also, the standard approaches cannot describe all the details of experiment and we hope that our algorithm can 
improve the situation. 

Let us finally note that in the 3D problems admitting full spherical symmetry (like spin glasses with random 
exchange interactions), the problem becomes reducible to ID case so that our algorithm might work in its present 
form. 

Acknowledgement: P. G. would like to thank Professor I. M. Sokolov for focusing our attention on the Gillespie's 
algorithm. All numerical simulations were completed by means of the facilities of the Platon - Science Services 
Platform of the Polish Pionier Network. 
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